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Abstract:  An  algorithm  is  presented  for  the  computation  of  the  second  funda¬ 
mental  tensor  V  of  a  Riemannian  submanifold  M  of  h3.  From  V  the  Riemann  curva¬ 
ture  tensor  of  M  is  easily  obtained.  Moreover,  V  has  a  close  relation  to  the 
second  derivative  of  certain  functionals  on  M  which,  in  turn,  provides  a 
powerful  new  tool  for  the  computational  determination  of  multiple  bifurcation 
directions.  Frequently,  in  applicatioru ,  the  manifold  M  is  defined  implicitly 
as  the  zero  set  of  a  submersion  F  on  Ir  .  In  this  case,  the  principal  cost  of 
algorithm  for  computing  V(p)  at  a  given  point,  p  ^  M  involves  only  the  decompo¬ 
sition  of  the  Jacobian  DF(p)  of  F  at  p  and  the  projection  of  d(d+l)  neighbor¬ 
ing  points  onto  M  by  means  of  a  local  iterative  process  using  DF(p) .  Several 
nvimerical  examples  are  given  which  show  the  efficiency  and  dependability  of 
the  method,  f  i 

J’  f 

J 

1 .  Introduction 


In  recent  years  various  computational  methods  for  the  analysis  of  dif¬ 
ferentiable  manifolds  have  been  developed,  see  e.g.  [8], [9], [10]  where  also 
other  references  can  be  found.  These  methods  have  numerous  applications  in  the 
study  of  multi -parameter  equilibrium  problems  and  their  bifurcation  behavior. 
But,  up  to  now,  there  appear  to  exist  no  general  purpose  numerical  methods  for 
the  computation  of  the  curvature  tensor  of  a  Riemannian  manifold  M  or  any  of 
its  related  quantities.  Yet,  the  curvature  tensor  incorporates  all  the  infor- 
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matlon  about  the  metric  properties  of  M  and  hence  is  of  considerable  impor¬ 
tance  in  the  solution  of  problems  which  intrinsically  involve  the  metric  of  M. 

In  this  paper  we  present  a  first  algorithm  for  the  computation  of  the 
second  fundamental  tensor  of  a  sub -manifold  M  of  r”  on  which  a  Riemannian 
metric  is  induced  by  the  Euclidean  inner  product  on  r”.  Our  reason  for  concen¬ 
trating  on  the  second  fundamental  tensor  is  two-fold.  First  of  all,  this  ten¬ 
sor  retains  all  the  information  about  the  metric  properties  of  the  manifold, 
and  the  Riemann  curvature  tensor  can  be  computed  from  it  by  means  of  a  simple 
formula  (see  e.g.[12]  and  section  2  below).  On  the  other  hand,  the  second  fun¬ 
damental  tensor  has  a  close  relationship  with  the  second  derivative  of  certain 
real  valued  functionals  on  the  manifold  and  this,  in  turn,  provides  a  powerful 
new  tool  for  a  computational  analysis  of  non -degenerate  bifurcation  phenomena 
when  possible  multiple  branching  occurs. 

In  Section  2  below  we  summarize  our  notation  and  collect,  without  proof, 
some  key  results  from  Riemannian  geometry.  Full  proofs  can  be  found  in  most 
texts  in  this  area,  as,  for  instance,  [1]  or  [12].  Section  3  Introduces  the 
concepts  leading  to  the  calculation  of  the  second  fundamental  tensor,  and  then 
in  Section  4  these  concepts  are  used  to  formulate  our  principal  algorithm.  In 
Section  5  we  discuss  the  indicated  application  of  this  algorithm  to  the  compu¬ 
tation  of  bifurcation  directions.  The  algorithm  has  been  implemented  as  a  gen¬ 
eral  purpose  FORTRAN  package  which  provides  also  an  interface  to  the  methods 
discussed  in  the  mentioned  earlier  references  [8] -[10]  as  well  as  to  the  con¬ 
tinuation  code  PITCON  (see  [11]).  Some  experimental  results  with  this  package 
involving  various  bifurcation  problems  are  presented  in  Section  6. 
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2.  Basic  concepts  from  Riemannian  geometry . 

We  denote  by  V  the  standard  connection  (covariant  differentiation)  in 
In  other  words,  for  any  smooth  vector  fields  X  and  Y  on  r’^,  is  the  vector 
field  obtained  by  differentiating  Y  in  the  direction  X. 

Throughout  this  paper,  M  C  r’^  will  be  a  smooth  sub -manifold  with  dimen¬ 
sion  d  >  2.  Then,  for  two  vector  fields  X  and  Y  that  are  tangent  to  M,  V^^Y  is 
defined  for  any  extension  of  X  and  Y  to  vector  fields  on  r’^.  As  a  vector  field 
on  M  (with  values  in  r”) ,  V^Y  turns  out  to  be  Independent  of  these  extensions. 

It  will  be  assumed  that  r"  is  equipped  with  the  canonical  inner  product 
<• , •>  which  induces  the  Riemannian  metric  on  M.  Accordingly,  orthogonality 
will  always  be  understood  in  the  sense  of  this  inner  product.  Then,  for  any 
vector  fields  X  and  Y  tangent  to  M  we  have 


7jjY  -  Vj^Y  +  V(X,Y)  , 


(2.1) 


where  V^^Y  and  V(X,Y)  are  the  tangential  and  normal  components  of  7^Y,  respec¬ 
tively.  In  other  words,  for  p  e  M,  (V.^Y)  and  (V(X,Y))  are  the  orthogonal 

A  p  p 

projections  of  onto  the  tangent  space  T^M  and  the  normal  space  N^M  of  M 

at  p,  respectively. 


While  the  operator  V  of  (2.1)  is  a  connection  on  M,  it  is  well  known 
that  V  is  a  symmetric  vector-valued  2 -covariant  tensor,  called  the  second  fun¬ 
damental  tensor.  In  particular,  the  value  (V(X,Y))p  depends  only  on  the 

values  X  and  Y  of  the  vector  fields  X  and  Y  at  p,  and  not  on  the  field 
P  P 

nature  of  X  and  Y. 


As  noted  in  the  introduction,  we  are  interested  in  the  numerical  evalua¬ 
tion  of  the  second  fundamental  tensor  V  at  a  given  point  p  6  M;  that  is,  in 
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the  computation  of  V(Xp,Yp)  for  an  arbitrary  pair  of  vectors  and  that 

are  tangent  to  H  at  p. 


As  observed  before,  the  second  fundamental  tensor  is  closely  connected 
with  the  Riemann  curvature  tensor  R  on  M.  In  fact,  for  any  four  vector  fields 
X,Y,Z,W  tangent  to  M,  the  following  simple  formula  holds 


<R(X,Y)Z,W>  -  <V(X,W) ,V(Y,Z)>  -  <V(X,Z) ,V(Y,W)>  ,  (2.2) 
(see  [12,  Vol.  IV,  p.  47]  and  recall  that  the  curvature  tensor  of  r’^  is  zero). 
In  particular,  for  W  -  X  and  Z  -  Y,  it  follows  that 


<R(X,Y)Y,X>  -  <V(X,X).V(Y,Y)>  -  |V(X,Y)|^  , 


(2.3) 


where  I*]  denotes  the  Euclidian  norm.  The  relation  (2.3), 
for  the  computation  of  the  curvature  k(P)  of  any  plane  P  - 
since 


in  turn,  allows 

span  {X  ,Y  )  C  T  M 
P  P  P 


where 


k(P) 


<R(X  ,Y  )Y  ,X  > 

p*  p*  p 

A(X  ,Y  ) 

P  P 


(2.4) 


A(X^,Y^)  - 
P  P 


2  2  2 
IX  I  |Y  I  -  <X  ,Y  > 

'  P'  '  P  P  P 


(2.5) 


is  the  square  of  the  area  of  the  parallelogram  generated  by  X^  and  Y^.  It  is 

well  known  that  k(P)  depends  only  on  the  plane  P,  and  not  on  the  specific 

choice  of  X  and  Y  . 

P  P 


In  order  to  clarify  the  mentioned  close  relationship  between  the  second 
fundamental  tensor  and  the  Hessian  of  real-valued  functions,  we  recall  first 
that  for  any  field  Z  of  normal  vectors  on  M,  the  second  fundamental  form  of  M 
in  the  direction  Z  is  the  bilinear  form 


(X,Y)  ^  <Vj^Z,Y> 


(2.6) 


defined  for  arbitrary  vector  fields  X  and  Y  tangent  to  M.  Because  of 
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<Z,Y>  "0,  it  follows  that 

0  -  X«<Z,Y>  -  +  <Z,VjjY>  -  <7^Z,Y>  +  <V(X,Y),Z>  . 

and  hence  that 

<V^Z,Y>  -  -<V(X,Y) ,Z>  (2.7) 

Now,  with  fixed  p  e  M  and  Z  e  N  M.  consider  the  functional 

7  :  M  R,  7(q)  -  <q  -  p,  Z^>,  for  q  e  M.  (2.8) 

It  is  a  standard  and  elementary  result  that  d7(p)  -  0;  that  is,  that  p  is  a 

critical  point  of  7.  The  Hessian  of  7  at  p  is  then  the  bilinear  form 

H^(Xp,Yp)  -  (X.(Y.7))p  .  (2.9) 

on  TpM,  where  X  and  Y  are  any  vector  fields  tangent  to  M  which  coincide  at  p 

with  X  and  Y  ,  respectively.  Let  Z  be  a  field  of  vectors  normal  to  M  such 

P  P 

that  Z  coincides  with  Z  at  p.  Then,  at  p,  the  second  fundamental  form  of  M 

P 

in  the  direction  Z  is  exactly  the  negative  of  the  Hessian  of  7  at  p  (see  e.g. 
[1,  p. 198]). In  other  words,  (2.7)  Implies  that 

H  (X  ,Y  )  -  <V(X  ,Y  ),Z  >  ,  (2.10) 

7  P  P  P  P  P 

In  Section  5  this  relation  will  provide  the  basis  for  the  determination  of 
bifurcation  direction. 

3.  Properties  of  the  second  fundamental  tensor. 

Let  X  and  Y  be  arbitrary  vectors  of  T  M  and  suppose  that 
P  P  P 

W  -  oX  +  BY  with  some  real  numbers  a  0  and  fl  »*  0 .  The  bilinearity  and 
P  P  P 

symmetry  of  the  second  fundamental  tensor  then  Implies  that 
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-  25  l''<''p'''p>  -  *^''<*p'V  ■  '’^'''’'p^’^p’!  • 

Obseirve  that  V(X  ,Y  )  can  be  computed  for  arbitrary  vectors  X  and  Y  of  T  M, 
P  P  P  P  P 

if  the  d(d  +  l)/2  values  V(X^,Xj),  1  <  i  <  j  <  d.  are  known  for  some  basis 
^  of  the  tangent  space  T^M.  Moreover,  because  of  (3.1),  the  computa¬ 
tion  of  the  quantities  V(Xj^,Xj)  can  be  reduced  to  that  of  the  d(d  +  l)/2  terms 
V(X^,X^),  1  <  i  <  d  and  (say)  V(Xj_  +  X^.X^  +  X^.),  1  <  i  <  j  <  d. 


Our  algorithm  for  the  calculation  of  V(X  ,Y  ),  will  be  based  on  an  ele- 

P  P 

mentary  geometric  construction.  Suppose  that  jX^j  -  1  and  introduce  the  space 

n  -  span(Xp)  ®  NpM  c  r”  .  (3.2) 

Note  that  dim  11  -  n  —  d  +  1  and  that,  because  of  the  (trivial)  relation 

r"  -  T  M  +  n  , 

P 

the  affine  space  p  +  n  intersects  the  manifold  M  transversally  at  p.  Thus, 
locally  near  p,  (p  +  H)  n  M  is  a  curve  C  and 


TC-nnTM-  span(X  )  . 

P  P  P 

Since  jX^I  -  1,  we  can  construct,  locally  near  p,  a  field  X  of  unit  vec¬ 
tors  tangent  to  M  which  along  C  is  tangent  to  C  and  coincides  at  p  with  X^. 
Indeed,  we  first  define  X  along  C  through  an  arclength  parametrization  of  that 
curve  and  then  extend  X  to  a  neighborhood  of  C  in  M.  Since,  by  construction,  X 
is  unit  vector  field  on  C,  it  is  trivial  that  this  extension  may  be  chosen  so 
that  |X|  -  1. 


Because  C  c  p  +  n  and  X  is  tangent  to  C  along  C,  it  follows  that 

e  n  along  C 


and,  in  particular,  that 


7 


e  n  .  (3.3) 

Moreover,  |X|  -  1  implies  that 

0  -  X*1  -  X*  <X,X>  -  2<^j^,X>  , 

whence  i-s  normal  to  X^  and  together  with  (3.2)  and  (3.3)  we  obtain 

<  v>p  ^  V_- 

In  other  words ,  the  tangential  component  of  vanishes  at  p  so  that  in  view 
of  (2.1). 

(7^)p  -  V(Xp.Xp)  .  (3.4) 

On  the  other  hand,  because  |X|  ~  1  and  X  is  tangent  to  C  along  C,  we  have  by 
definition 


|(Vj^)p|  -  kp  (>  0)  ,  (3.5) 

where  is  the  curvature  of  C  at  p.  Moreover,  if  k^  0,  the  principal  nor¬ 
mal  n^  of  length  jn^l  -  1  is  defined  and 

<v^p  -  Vp  ’ 

Altogether,  (3.4),  (3.5)  and  (3.6)  show  that 


V(X^,X^)  -  0  if  k^  -  0  , 

P  P  P 

V(X  ,X  )  -  k  n  if  k  >  0  . 
P  P  P  P  P 


When  k  »<0 , 
P 


the  osculating  plane  to  C  at  p  is  the  plane 


span  {X  ,n  )  e  r"  . 

P  P 

A  standard  result  of  the  local  theory  of  curves  states  that  the  osculating 

plane  is  the  limit  of  the  planes  span(qj^-p,  q2~p)  for  any  two  distinct  points 

p,  i“l,2  of  C  that  tend  to  p.  Since  k^  0  it  follows  that  p  and  q^^,  q^ 

3 

are  not  co-linear  (see  [12,  Vol.  II]  for  curves  in  R  ;  the  argument  extends 
easily  to  curves  in  r”)  .  Thus  qj^,  q2  and  p  define  a  circle  in  r”  and  for 
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-►  p,  i-1,2,  the  limit  of  these  circles  is  a  circle  C  in  the  osculating 
plane  with  radius 


r  -  lA  (3.8) 

P  P 

This  is  the  osculating  circle  at  p . 


4.  A  Computational  Algorithm  for  the  Fundamental  Tensor 

In  general,  whenever  manifolds  arise  in  numerical  calculations,  they  are 
defined  either  through  a  local  parametrization  which  is  an  immersion  or, 
implicitly,  as  the  zero  set  of  a  submersion.  As  in  [8] -[10]  it  will  suffice  to 
consider  only  the  latter  case.  Hence,  suppose  that  F:  n  c  r”  -►  R™  is  smooth  on 
the  open  set  Q  and  rank  DF(q)  -  m  for  all  q  e  0.  Ue  assume  always  that 
d  -  n-m  >  2.  Then 


M  -  {q  e  0  ;  F(q)  -  0) 


(4.1) 


is  a  d*dlmenslonal  submanifold  of  r".  For  any  p  e  M  we  identify  the  tangent 
space  TpM  with  ker  DF(p)  and  the  normal  space  N  M  with 
(ker  DF(p))  -  rge  DF(p)  where  as  usual  the  asterisk  denotes  the  adjoint 
operator.  There  are  various  possible  methods  for  the  computation  of  the  ortho¬ 


normal  bases  of  T^H  and  for  Instance,  a  simple  approach  is  based  on  the 

use  of  the  QR- factorization  of  DF(p)  ,  (see  also  [10]). 


For  any  given  €  T^M,  ^  0,  and  sufficiently  small  nonzero  h  G  R  the 
implicit  function  theorem  ensures  the  existence  of  a  unique  point  w(h)  e  N^M 
such  that 


F(p  +  hXp  +  w(h))  -  0  . 

Once  again,  there  are  various  methods  for  computing  w(h) . 


(4.2) 

In  particular,  if 
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the  QR* factorization  of  DF(p)  Is  already  available,  a  chord-Gauss 'Newton  pro¬ 
cess  can  be  applied  --  with  q  -  p  +  hX^  as  starting  point  --to  project  q  onto 
M  along  Ue  refer  to  [10]  for  details  and  a  convergence  result.  Note  that 
the  process  of  projecting  q  onto  M  requires  only  the  evaluation  and  decomposi¬ 
tion  of  the  Jacobian  of  F  at  p.  Of  course,  we  may  use  other  decompositions 
than  the  QR- factorization. 


As  in  the  previous  section  and  with  the  notation  (3.2),  let 


C  -  M  0  (p  +  n)  denote  the  curve  induced  on  M  by  the  given  vector  X  €  T  M. 

P  P 

Then,  for  small  |h|  we  have  p  +  hX^  +  w(h)  e  M  n  (p  +  II)  -  C.  This  suggests 
that  we  compute  for  some  h  >  0  the  points 


qj^  -  p  +  hXp  +  w(h)  ,  q2  -  P  +  hX^  +  w(-h) 
on  M.  Then,  if  0  and  h  is  sufficiently  small,  the  triangle  formed  by 

q^^,  i-1,2,  and  p  in  the  plane  span(q^-p,  q2-p)  is  not  degenerate. 

The  well-known  Heron  formula  of  planar  geometry  states  that  for  a  non¬ 
degenerate  triangle  with  side  lengths  a,b,c,  the  curvature  k  of  the  cir¬ 
cumscribing  circle  is  given  by 

k  -  ^|^[s(s-a)(s-b)(s-c)]^'^^  (4.3) 

where  s  -  (a+b+c)/2.  In  our  case,  the  sides  a  -  [q^^-pl  and  b  -  |q2-p|  are 

2 

equal  up  to  order  h  and  the  third  side  is  almost  equal  to  their  sum.  Hence  it 
is  natural  to  introduce  the  scaled  quantities  a^  -  a/c ,  and  b^  -  b/c,  and  to 
rewrite  (4.3)  in  the  form 

-111  2  h  2  h 

k  -  ip  +  ^]  (1  -  «  )  (1  -  7  )  (4.4) 

c  c 

where  5  -  a^  -  b^  and  y  -  l/(a^  +  b^).  The  term  in  square  brackets  is  close  to 
4  while  the  first  square  root  is  approximately  equal  to  one.  Thus,  both  of 
these  terms  can  be  evaluated  safely.  But,  in  the  last  square  root  y  is  close 


/ 
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to  one,  and  hence  we  compute  that  term  as 

^  -  arcos(7),  (1  -  7^)^  -  sin(^)  (4.5) 
When  1-7  falls  below  machine  precision  then,  in  floating  point  arithmetic,  tj) 
will  be  zero  and  we  set  k  -  0. 


Clearly,  when  h  tends  to  zero  then  the  circle  through  the  points 
i-1,2  and  p  tends  to  the  osculating  circle  and  hence  k  becomes  the  curva¬ 
ture  kp  of  the  curve  C  at  p.  Thus  our  algorithm  produces  an  approximation  k  of 

k  . 

P 

If  k  is  not  zero,  then  an  approximation  n  of  the  normal  vector  n^  can  be 
generated  by  orthogonalizing  the  vector  v  -  (q^  -  p)  +  (q2  “  P)  with  respect 
to  Xp  and  then  normalizing  the  result  to  length  one;  that  is,  by  applying  the 
algorithm 


Thus  altogether, 

V(X  X  ). 

P  P 


we 


n  V  -  <X^,  v>  X  ,  n  : - — 

^  P  |nl 

have  obtained  an  approximation  V(X  ,  X  )  -  k  n 

P  P 


of 


As  noted  in  the  previous  section,  the  calculation  of  the  fundamental  ten¬ 
sor  V  can  be  reduced  to  the  computation  of  the  d(d  +  l)/2  quantities 
V(X^,X^),  1  <  i  <  d  and  V(,X^  +  Xj,Xj^  +  Xj),  l<i<j<d  for  some  basis 
Xj^,  i  —  1,  .  .  .  ,d  of  TpM.  Of  course,  this  basis  may  be  chosen  in  various  ways; 
we  indicate  here  only  some  advantageous  choices  in  the  cases  d  -  2,3.  For  this 
purpose,  suppose  that  an  orthogonal  basis  of  T^M  is  already  available  which 
then  defines  an  isomorphism  U  from  onto  TpM. 

2 

For  d  -  2  we  introduce  in  R  the  vectors 
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4 

where,  for  ease  of  notation,  r  -  3  .  Then, 
values  -  V(Xj^,  ,  i  -  1,2,3,  of  the 
computed,  and,  by  (3.1),  we  obtain 


^3- 


X1  +  X2 


with  the  above  algorithm,  the 
second  fundamental  tensor  can  be 


V12  -  V(X^,  X2)  -  +  ^2  -  V3)  . 

Now  for  any  two  vectors  ^2  €  T^M  it  follows  that 

Yf  -  .  i  -  1.2 

with 


“1 

2 

’2  1' 

"  3 

1  2 

CM 

and  hence  that 


V(Yi.  Y2)  -  0^02  +  (a^^2  ^  “2^1>'^12  ^  Wz 

For  d  -  3  we  proceed  analogously  and  introduce  the  vectors 


X^  -  u 

1' 

0 

.  Xj  -  1  u 

-1' 

T 

■  *3  -  4  ” 

-1' 

— r 

.0. 

.  0. 

.2r 

and 


(4.6) 


X^  -  X^  +  X2  .  X5  -  +  X3)  ,  Xg  -  (|)^(X2  +  X3) 

Then,  once  again,  we  compute  -  V(Xj^,  X^^) ,  i  -  1,...,6,  and 

V(X^,  X2)  -  V^2  “  "  ’^l  "  ^^2^  ’ 

''«!•  Jts)  -  ''13  -  I’s  -  2<''l  ''3>  • 

VCXj,  X3)  -  V23  -  fVj  -  icVj  ..  V3)  . 
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Now  for  any  three  vectors  Yj^eTM,  i-  1,2,3,  we 


have 


with 


and  hence 


Yi  -  +  /9^X2  +  7j^X3  ,  i  -  1,2,3 


f  '\ 

T 

YtX- 

“i 

'5  3  2' 

i  1 

^i 

1 

“  3 

3  5  2 

T 

V2 

2  2  4 

T'i 

^  J 

yTx- 

V  J 

i  3 

V.  J 

V(Y^,  Y^)  -  "^1^3 


.  i  -  1.2,3 


+  “2^1^^12  ‘‘‘  ^^1^2  2Z  ■*■  ^“l'’^2  “2”’'l^^l3' 


It  should  be  evident  how  we  might  proceed  for  higher  dimensional  mani¬ 
folds  M.  Note  that,  independent  of  the  dimension  d  of  M,  the  computation  of 
the  fundamental  tensor  at  a  point  p  of  M  requires  only  one  evaluation  and 

decomposition  of  the  Jacobian  of  F,  namely  at  p.  The  computational  cost  for 

3 

this  is  of  order  n  ,  where,  of  course,  n  is  the  dimension  of  the  embedding 

space.  Once  the  decomposed  Jacobian  at  p  is  available,  the  cost  of  projecting 

2 

each  one  of  the  d(d+l)  required  points  onto  M  is  of  order  n  while  all  other 
parts  of  the  algorithms  involve  a  lower  order  of  operations. 


5.  De termination  of  Bifurcation  Directions 

As  in  section  4,  suppose  again  that  the  smooth  mapping  F:  0  c  -♦  R™, 
(d  ••  n-m  >  2)  satisfies  rank  OF<q)  -  m  for  all  q  on  the  open  set  0  and  hence 
that  (4.1)  defines  a  d-dlmenslonal  submanifold  of  r".  At  a  given  point  p  €  M, 
we  consider,  as  in  [8],  local  coordinate  systems  which  are  induced  by  some  d- 
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dimensional  linear  subspace  T  of  R^.  More  specifically,  if  T  n  N^M  -  (0)  then 
it  follows  from  the  Implicit  function  theorem  that  there  exist  open  neighbor¬ 
hoods  Q  c  T  of  the  origin  of  T,  and  S  c  r”  of  p,  respectively,  and  a  smooth 

mapping  f;  Q  -»  T"*"  such  that  f(0)  -  0  and  M  n  S  -  0(Q)  where 

0:  Q  -»  r",  0(u)  -  p  +  u  +  f(u)  e  r",  u  e  Q  .  (5.1) 

Note  that  always  D0(O)u  e  T  M  for  u  €  Q  and  hence  that  for  the  local  coordi- 

P 

nate  system  induced  by  T  -  T^M  we  have  Df(0)  -  0.  For  this  case  we  consider 
here  the  numerical  evaluation  of  the  bilinear  mapping 


(X  ,Y  )  €  T  M  X  T  M 
P  P  P  P 


D^f(0)(X^,Y  )  . 
P  P 


(5.2) 


Let  Z.,  ,  ...  ,Z  be  an  orthonormal  basis  of  the  normal  space  N  M  and  hence 

1  -  m  p 

f  -  f^  Z^  +  ...  +f^  Z^,  f^  -  <f,  Z^>,  i^l . m.  (5.3) 

For  each  Z^,  1  <  i  <  m,  we  consider,  as  in  (2.8),  the  functional 

M  R,  7j^(q)  -  <q  -  p,  Z^>.  for  q  £  M  (5.4) 

for  which 


7j^(q)  -  <u  +  f(u),Zj^>  -  <f(u),Zj^>  -  fj^(u),  for  q  e  Q,  (5.5) 

and  hence  f  -  y.*9.  Since,  Df(0)  -  0  we  have  D7. (0)  -  0  and  D0(O)X  -  X  for 
11  1  P  P 

any  X^  e  T^M.  Thus  the  Hessian  of  7^^  at  p  is  given  by 

H  (X  Y  )  -  D^f  (0)(X  Y^),  for  all  X^,Y^  e  T  M  . 

Ipp  Ipp  PPP 

Now  (2.10)  implies  that 


D^f  (0)(X^,Y  ) 

i  p  p 

whence,  by  (5.3), 


<V(X^,Y^),Z.>,  i-1, 

p  p  i 


,m. 


(5.6) 


D^f(0)(X  Y^)  -  V(X„,Y^). 

P  P  P  P 

In  other  words, the  second  fundamental  tensor  of  M  at  p  is  exactly  the  second 
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derivative  of  f  at  the  origin  of  T^M.  Note  that  this  simple  relation  is  valid 
only  when  Df(0)  -  0;  that  is,  when  the  local  coordinate  system  at  p  is  induced 
by  the  tangent  space.  But  that  is  a  natural  choice  in  the  framework  of  bifur¬ 
cation  problems. 

Generally,  in  applications  a  d-dimensional  natural  parameter  subspace 
A  c  r"  is  identified.  We  follow  here  the  approach  in  [4]  and  call  a  point 
p  €  M  a  foldpoint  (with  respect  to  A  )  if  A  does  not  induce  a  local  coordinate 
system  at  p;  that  is,  if 

Nrt  -  A  n  N  M  »*  (0) . 

0  p 

The  integer  r^^  -  dimCN^)  is  the  first  singularity  index  of  the  foldpoint  and 
the  set 

(q  €  M;  <q  -  p,  n>  -  0,  for  all  n  e  Nq) 
the  cutset  of  H  at  p.  Let  Z^,  i-l,...,rj^  be  an  orthonormal  basis  of  Nq  and 
extend  it  to  an  orthonormal  basis  of  all  of  N^M.  Then  it  follows  from  (5.4) 
and  (5.5)  that,  locally  near  p,  the  cutset  is  the  zero  set  of  the  mapping 


g:  Q  c  TpM 


R  .  g(u)  -  (fj^(u), 


,f  (u))' 


u  e  Q. 


(5.7) 


Evidently,  we  have  g(0)  -  0,  Dg(0)  -  0,  and  if  r2  >  1  exists  such  that 


D‘"g(0)  -  0,  k-0,1 . r2,  D  g(0)  0 

then  r2  is  called  the  second  singularity  index  of  p.  Under  some  non- degeneracy 
conditions  (see  e.g.  [2], [4],  or  [6]),  the  form  of  the  cutset  is  determined, 
locally  near  p,  by  the  nontrivial  zeroes  of  the  (r2+l)-form 

tz+l 

g(0)(X^ . X^),  X^  6  T  M  . 

P  P  P  P 

For  r2  -  1  it  follows  from  (5.7)  that 
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D‘g(0)(X„.X^)  -  (<V(X^.X^).Z,> . <V(X„,X^).Z^  »,  X^  6  T  M.  (5.8) 

PP  PP'*-  PP  P  P 


P  P 

In  other  words ,  the  bifurcation  directions  at  p 


are  the  zeroes 


€  T  H,  X  0  of  the  r^  quadratic  equations 


<V(Xp.Xp).Zj^>  -  0.  i-1 . r^. 


(5.9) 


It  may  be  useful  to  reformulate  this  result  for  the  standard  setting  of 
bifurcation  theory.  For  this  suppose  that  G:  R™xR^  R™  is  a  smooth  mapping  on 
some  open  neighborhood  of  a  point  (y,A)  e  r’^'^xR^  where  G(y,A)  -  0  and 
dim  ker  DG(y,A)  -  r+1  >  2.  For  ease  of  notation  we  shall  write  here  also 
X  -  (y,A)  and,  in  particular,  x  -  (y,A).  Note  that  the  conditions  on  G  imply 

that  dim  ker  DyG(x)  -  r  >  1  and  dim  rge  DG(x)  -  m-r  <  m-1.  Let  a^^ . a^  e  R™ 

be  linearly  independent  vectors  which  span  ker  DG(x)  .  Then  span(a^ . a^)  is 

a  complement  of  rge  DG(x)  in  R™  which  suggests  the  introduction  of  the  unfold¬ 
ing 


F:  R®xA  -♦  R*^,  A >  R^xR®,  F(y,A,fi)  -  G(y,A)  +  J^^a^  +.  . .+  fi^a^  (5.10) 
where  S  -  (i , 5^) .  With  n-m+l+r,  d-r+1,  q  -  (y,A,fi)  this  corresponds 

e 

exactly  to  the  earlier  setting.  Clearly,  we  have  F(p)  -  0  at  p  -  (x,0)  and 
the  condition  rge  DF(p)  -  m  holds  in  an  open  neighborhood  of  p.  Hence  the  d- 
dimensional  manifold  (4.1)  is  well-defined  and  p  is  a  foldpoint  with  respect 
to  A  with  first  singularity  index  r^^  -  r.  Since 


DF(p)  I-l . r 

where  ej  are  the  natural  basis  vectors  of  r”,  the  set  N^  -  A  nN^M  is  spanned 
here  by  . ^  ®nd  the  cutset  is  the  solution  set  of  the  equa¬ 
tions  F  -  0,  S  -  0  and  hence  of  G  -  0.  Our  local  coordinate  system  is  given 
now  by 


e(e)  -  (x+e,5(e)),  4  e  ker  DG(x) 
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Thus  we  have  t-l,...,r  and  (5,6)  becomes 

D^«i(0)(|.»j)  -  <V((e.0),(»,.0)),Z^>.  i-1 . r,  for  ^,.7  e  ker  DG(x) . 

Hence,  If  Che  foldpolnt  has  second  singularity  index  r^  -  1,  then  under  Che 
mentioned  non- degeneracy  condition  the  bifurcation  directions  for  G  -  0  are 
exactly  the  solutions  of  the  quadratic  equation 

D^6(0)(|.^)  -  0  (5.11) 

In  bifurcation  theory,  the  characterization  of  the  bifurcation  direc¬ 
tions  ,  given  here  in  terms  of  6 ,  is  usually  formulated  in  terms  of  the  reduced 
mapping  obtained  after  transforming  the  equation  G  »  0  by  means  of  the 
Lyapunov- Schmidt  reduction.  While  both  characterizations  give  equivalent 
results,  Che  method  we  use  here  completely  bypasses  Che  Lyapunov- Schmidt 
reduction  and  allows  for  all  calculations  to  be  performed  in  the  framework  of 
the  smooch  manifold  (4,1)  on  which  no  singularity  hampers  the  computations. 

The  equations  (5.9)  (or  (5.11))  determining  the  bifurcation  directions 
are  readily  solvable.  In  fact,  suppose,  as  in  section  4,  that  X.,  , . .  .  ,X .  is  a 
basis  of  TpM  for  which  the  components  -  V(Xj^,Xj),  l<i<j<d  of  the 

second  fundamental  tensor  have  been  computed.  Then  (5.9)  reduces  to  a  system 
of  homogeneous  quadratic  equations 

x^A^x  -  0,  x€R,x>‘0  (5.12) 

where  the  dxd-matrices  have  the  elements 

"  ^iJ’V’  1  <  i.j  <  d.  k  -  1 . 

and 


X  -  x-X,  + 

p  11 


This  system  of  quadratic  equations  can  have  finitely  many  isolated  solutions 


17 


only  in  the  case  d  -  r^-t-1.  Hence  this  is  the  only  case  for  which  solutions 
will  be  sought  numerically.  Clearly,  for  small  d  the  system  is  easily  solv¬ 
able.  The  case  d  -  2,  (r^  -  1}  is,  of  course  trivial,  and  also  for  the  case 
d  -  3,  (r^  -  2)  simple  methods  are  available.  In  fact,  (5.12)  has  non-trivial 
solutions  only  when  all  matrices  are  indefinite.  Thus  in  the  (3x3) -case  each 
must  have  one  eigenvalue  of  different  sign  than  the  other  two.  Hence,  if 
there  is  no  degeneracy,  then,  after  diagonalization,  the  first  of  the  two 
equations  has  the  fonn 

^1^1  “  ^2^2  ^*3^3 

with  fi^>  0,  i-1,2,3.  This  suggests  the  normalization  -  1  and  the  use  of 
the  elliptic  coordinates  $2  ”  “  cos(t),  ^  sin(r),  0  <  r  <  rr  where 

a  -  ,  P  -  (/i^Z/ij)  .  Now  we  can  readily  detect  the  intervals  in  the 

r -variable  on  which  the  corresponding  residual  of  the  second  equation  changes 
sign,  which,  in  turn,  allows  the  application  of  a  standard  root  finder  for 
determining  the  solutions. 

Clearly,  for  larger  dimensions  such  simple  approaches  are  not  readily 
available,  and  the  system  (5.12)  has  to  be  solved  by  means  of  one  of  the  poly¬ 
nomial  root  finders  as,  for  example,  the  CONSOL  system,  [5j. 


6.  Ntunerical  Examples 

In  order  to  demonstrate  the  performance  of  the  algorithm,  we  present  here 
some  numerical  results  obtained  with  a  FORTRAN  implementation  run  in  double 
precision  on  a  VAX-cluster. 

As  a  first,  very  simple  problem,  consider  the  cusp 
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F:  -*  R^,  F(x)  -  -  X2.*2  ”  *3’  x  €  R^  .  (6.1) 

The  tangent  space  at  x  -  0  is  T^M  -  span(ej^,e2)  where  e^^,  i-1,2,3  are  the 

3 

natural  basis  vectors  of  R  .  A  straightforward  calculation  shows  that  at  x  -  0 
we  have 

V(u^,U2)  -  -sih(^j^  +  ^2)®3'  “  (cos(^j^)  ,sin(^j^)  ,0)^  e  T^M,  i-l,q5.2) 

Table  1  gives  some  comparison  of  computed  values  of  the  curvature  and  their 
corresponding  exact  values  -  |sin(^^  +  ^2) I  different  choices  of  the 

two  angles.  Here  we  used  a  step  h  -  0.005  in  the  curvature  evaluations  and  a 
tolerance  of  10  ^  for  the  projection  of  the  points  onto  the  manifold. 


1  Table 

1 

1:  Curvature  of  the 

Cusp  Function  | 

1 

1 

Exact 

Computed  j 

1 

1 0 

1 

0 

0.0 

0.0  1 

I 

1 0 

1 

0.3x 

0.809017 

0.809013  1 

1 

1 0 

0 . 6ir 

0. ‘951057 

0.951052  1 

1 

1 0 

1 

0.9ir 

0,309017 

0.309016  1 

1  0.3jr 

1 

0.3ir 

0.951057 

0.951052  1 

1  0.3x 

1 

0.4ir 

0.809017 

0.809013  1 

1  0.3w 

0.5ir 

0.587785 

0.587782  1 

1 

1  0.6jr 

1 

0.6ir 

0.587785 

0.587782  1 

1 

1  0 . 6jr 

0.9x 

1.0 

0.999995  1 

1 

1  0.8x 

j 

0.9ir 

0.809017 

0.809013  1 

The  computed  values  certainly  are  in  excellent  agreement  with  the  exact  data. 
All  our  experience  so  far  Indicates  that  the  curvature  algorithm  is  Indeed 
very  reliable,  not  just  for  small  problems  as  this  one. 
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As  a  second  example,  consider  the  two  point  boundary  value  problem 

-u"  -  Xu  +  au^  -  0,  u(0)  -  u(jr)  -  0  (6.3) 

which,  after  a  straightforward  discretization  and  In  unfolded  form,  leads  to 
the  finite  dimensional  system 

Ay  +  h^(aQ(y)  -  Ay)  +  /iw  -  0,  y  6R^,  h  -  >r/(k+l)  .  (6.4) 

2  2  T 

Here  Q(y)  -  (y^^ . y^^)  and  A  Is  the  symmetric,  trl-dlagonal ,  kxk  matrix  with 

diagonal  elements  2  and  sub-diagonal  entries  -1.  Thus,  with  q  -  (y,A,/x)  € 
and  n  -  k+2,  m  -  k,  d  -  2,  (6.4)  Is  a  problem  of  the  form  considered  In  sec¬ 
tions  4  and  5 . 


Let  <^2^  <  •  •  •  <  be  the  eigenvalues  of  A.  The  unfolding  vector  w  was 

chosen  as  a  normalized  eigenvector  of  A  corresponding  to  the  eigenvalue  for 

2 

a  given  1.  Then  p  -  (0,a^/h  ,0)  Is  a  foldpolnt  of  (6.4)  with  respect  to  the 

parameter  space  spanned  by  the  basis  vectors  e  - ,  e  of  and  both  slngular- 

n— i  n 

Ity  Indices  of  p  equal  one.  A  theoretical  analysis  of  (6.3)  shows  that  each  p 
Is  a  bifurcation  point  and  that  for  odd  values  of  1  the  bifurcation  Is  tran- 
scrltlcal  while ‘for  even  1  the  branches  Intersect  at  a  right  angle.  More 
specifically,  for  odd  1  the  angle  between  the  bifurcating  branches  satis¬ 
fies 


-  cos(aj^)  -  r/(l  +  r^)^, 
where  w^^,  .  .  .  ,w^  are  the  components  of  w. 


(6.5) 


Table  2  shows  the  computed  and  theoretical  values  of  for  1  -  l,...,k 
In  the  case  k  -  15.  Again  the  agreement  between  the  computed  and  the  predicted 
data  Is  excellent. 
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Table  2:  7 -values  (6.5) 

1  Exact  Computed 


1 - 

1 1 

0.287444 

0.287445  1 

1  2 

0.0 

0.155(-15)  1 

1  3 

0.996991 

0.996700  1 

1  ^ 

0.0 

0.673(-16)  1 

!  5 

0.608096(-l) 

0.608102(-1)  1 

1  6 

0.0 

0.0  1 

1  7 

0.462437(-l) 

0. 462441 (-1)  1 

1  8 

0.0 

0.171(-11)  1 

1  9 

0.478171(-1) 

0.478175(-1)  1 

1  10 

0.0 

0.160(-15)  1 

1  11 

0. 940416 (-1) 

0.940424(-l)  1 

1  12 

0.0 

0.801(-17)  1 

1  13 

0.340802(-2) 

0.340805(-2)  1 

1  1^ 

0.0 

0.171(-11)  1 

1  15 

1 - 

0.869785(-4) 

0.869794(-4)  | 

.  1 

As  a  third  problem  we  consider  the  discrete  Brusselator  (see  e.g.  [7])  in 
unfolded  form 

Ay  -  Ah^((^-l)y  +  9z  +  G(y,z,^))  +  -  0 

2  h  -  jr/(k+l),  y,z  e  (6.5) 

Ay  +  0.25Ah  ()8y  +  9z  +  G(y,z,/3))  +  i/w^  -  0 

where  again  A  is  the  tri-diagonal,  kxk  matrix  used  in  (6.4),  and 


G(y,z,^)  -  . g(yj^,z^j,^))’^,  g(s,t,/3)-(|^  +  t)s^  +  6st 


With  q  -  (y,z, A, /!,»/)  e  R*'  and  n  -  2k+3,  m  -  2k,  d  -  3  this  is  a  problem  of  the 
form  considered  in  sections  4  and  5.  If  again  <.  .  .  <  sctb  the  eigenvalues 
of  A  then  some  calculation  shows  that  p  —  (0,0, A, 0,0)  with 


A  -  +  i  j  (6-6) 
3h  J 

is  a  foldpoint  of  (6.5)  with  respect  to  the  parameters  A,  /i,  v  and  with  first 
and  second  singularity  indices  2  and  1,  respectively.  For  k  -  3  and 
i  -  1,  j  -  2,  (6.6)  gives  the  values  A  -  1.08158,  /3  -  6.83343  used  in  [7], 
while  for  k  -  8  they  become  A  -  1.14009,  ^  -  6.96599.  In  each  case,  our  algo¬ 
rithm  determined  four  bifurcation  directions.  For  k  -  3  the  computed  direc¬ 
tions  are  given  in  Table  3.  This  appears  to  be  the  first  computational  deter¬ 
mination  of  bifurcation  directions  for  singularities  with  co-dimension  larger 
than  2  in  the  literature. 


Table  3:  Bifurcation  Directions  for  the  Brusselator 
#1  #2  #3  #4 


!  -0.646189(-6) 

0.777229 

0.253517 

-0.380313 

-0.594689(-6) 

0.370361 

-0.370326 

-0.537862 

-0.l94827(-6) 

-0.253460 

-0.777237 

-0.380338 

-0.311403(-6) 

-0.321441 

-0.292234(-l) 

0.212202 

-0.331812(-6) 

-0.206646 

0.206627 

0.300105 

-0.157850(-6) 

0.291988(-1) 

0.321438 

0.212211 

-1.000000 

0.218236 

-0.218217 

-0.491205 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

Clearly,  once  the  bifurcation  directions  are  available,  a  standard  continua¬ 
tion  process  (see  e.g.  [11])  can  be  started  to  follow  the  corresponding  solu¬ 
tion  branch. 

The  examples  certainly  Indicate  that  th^^lgorithms  presented  here  per¬ 


form  very  well  in  practice.  As  noted,  our  experience  has  shown  that  the  curva¬ 
ture  computations  are  very  reliable.  Of  course,  in  the  above  examples  we 
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always  worked  with  an  exactly  known  foldpolnt  p  with  known  singularity 
Indices.  In  practice,  one  often  has  only  some  approximation  of  such  a  point. 
But  in  our  approach,  this  is  not  a  serious  concern.  In  fact,  the  second  funda¬ 
mental  tensor  V  is  a  smooth  function  of  the  base  point  on  the  manifold  H  and 
hence  the  computed  V  can  be  expected  to  be  an  approximation  of  V  at  the  exact 
foldpolnt.  Thus  for  any  chosen  normal  vectors  the  quadratic  equations  (5.9) 
should  also  be  close  to  the  corresponding  equations  at  p  itself.  This  leaves 
only  the  determination  of  the  critical  first  singularity  index  r^^  and  of  the 
normal  vectors  that  span  the  intersection  A  n  N^M.  We  refer  to  [7]  for 
some  comments  about  the  determination  of  r^,  and  to  [3]  for  a  method  to  com¬ 
pute  an  approximation  of  the  basis  vector  in  the  case  r^  -  1. 
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